Tutorial 9: Advanced Modelling and Visualisation

Programming for Engineers (ENGG1001)

Published

February 18, 2026

Tutorial Tasks

Task 1: Creating a Meshgrid for 3D Surface Plots

So far, you have had practice with creating 1D numpy arrays and plotting them. For 3D surface plots, we need to create a grid of points in the x-y plane. Write a function create_meshgrid(length: float, num_points: int) -> tuple[np.ndarray, np.ndarray] that creates a meshgrid of points from -length to length with num_points for both x and y axes. For a further explanation of what a meshgrid is, see the callout below.

What exactly is a meshgrid? A meshgrid is a 2D grid of points that can be used for plotting surfaces or evaluating functions over a 2D domain. It consists of two 2D arrays. Consider the example below where we create a 10 x 10 meshgrid:

x_mesh:
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
 [0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]
y_mesh:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
 [3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
 [4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
 [5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
 [6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
 [7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
 [8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
 [9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]

Essentailly, a meshgrid is a way of presenting a set of cartesian coordinates in a structured format. So if we feed an index (which essentially acts as a coordinate) into the meshgrid, we can get the corresponding x and y values. For example, if we want to get the x and y values at the index (2, 3), we can do:

x value at x_mesh[2, 3]: 3.0
y value at y_mesh[2, 3]: 2.0

Here is a sample usage of the function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
print(x_mesh.shape)
 
(100, 100)
>>>
print(x_mesh[0, 0])
 
-10.0
>>>
print(y_mesh[50, 0])
 
3.303030303030301
Tip

What function could you use to create the 1D arrays for x and y? Once you have created the 1D arrays, then you can use numpy.meshgrid to create the 2D grid.

Task 2: Creating a Plane Surface

Now that we can create a meshgrid of points, we can use it to create a 3D surface. Write a function plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray that computes the z values for a plane surface defined by the equation: \[z = ax + by + c\]

Here is a sample usage of the function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)
>>>
print(plane.shape)
 
(100, 100)
>>>
print(plane[0, 0])
 
-27.0
>>>
print(plane[50, 50])
 
3.303030303030301
Tip

To compute the wave surface, you will need to calculate the radial distance R for each point in the meshgrid and then apply the RBF formula. Remember to use vectorisation!

Task 3: Modelling waves with a Radial Basis Function

Now that we have a meshgrid of points, we can model a wave surface using a Radial Basis Function (RBF). The RBF is defined as: \[f(x, y) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r}\] Where:

  • \(A\) is the amplitude of the wave
  • \(k\) is the wave number
  • \(r\) is the radial distance from the origin: \(r = \sqrt{x^2 + y^2}\)
  • \(\omega\) is the angular frequency
  • \(\alpha\) is the decay constant of the wave with respect to distance
  • \(t\) is time in seconds

Write a function radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray that computes the wave surface values for the given parameters. Here is a sample usage of the function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
>>>
print(wave_surface.shape)
 
(100, 100)
>>>
print(wave_surface[50, 50])
 
0.9458561801297212
>>>
print(wave_surface[0, 0])
 
-0.24310473049518086

Task 4: Visualising 3D Surfaces

Now that we have a function to compute the wave surface, we can visualise it using Matplotlib’s 3D plotting capabilities. Write a function plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray) that creates a 3D surface plot of the surface.

Here is a sample usage of the function potting the plane surface:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)
>>>
plot_surface(x_mesh, y_mesh, plane)

Which produces the following 3D surface plot:

3D Surface Plot of the Plane Surface

Also try plotting the wave surface using the same function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
>>>
plot_surface(x_mesh, y_mesh, wave_surface)

Which produces the following 3D surface plot:

3D Surface Plot of the Radial Wave
Tip

To be able to setup an axis for a 3D plot, you need to specify the projection argument when creating the subplot. The plot_surface method can be used to create the surface plot. You can use the code below to get started:

fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.plot_surface(X, Y, Z)

If you want to make your plot look nicer, you can also specify a colormap using the cmap argument in plot_surface. For example, you can use cmap='viridis' to use the Viridis colormap, like the example above.

Note

Now that you have created an easy way to visualise the wave surface, determine what effects the following parameters have on the wave surface by changing their values and observing the resulting plot:

  1. Amplitude A
  2. Wave number k
  3. Angular frequency omega
  4. Decay constant alpha
  5. Num points in the meshgrid num_points

Task 5: identifying local maxima and minima

Now that we have a visualisation of the wave surface, we can identify the local maxima and minima of the surface. Write a function find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray] that finds and returns the radius values (how far from the origin) of the local maxima and minima of the surface. Here is a sample usage of the function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
>>>
maxima, minima = find_local_extrema(wave_surface)
>>>
print(maxima)
 
[9.39448244 6.16244406 3.1329419 0.14284985 3.1329419 6.16244406
 
9.39448244]
>>>
print(minima)
 
[7.77843366 4.7485492 1.51851479 1.51851479 4.7485492 7.77843366]
Tip

To find the local maxima and minima of radial function, we can use the fact that the surface is “radially symmetric”, meaning that the surface values only depend on the radius from the origin. As such, we can compute the radius values for a central “slice” of the surface (e.g. a given row or column of the surface), and then find the local maxima and minima of that slice.

To illustrate this, consider the following central slice of the curve above:

Central Slice of the Radial Wave Surface

We can see that the local maxima and minima of the surface correspond to the local maxima and minima of this slice. Therefore, we can find the local maxima and minima of this slice to determine the radius values of the local maxima and minima of the surface.

To find the local maxima and minima of a 1D array, you can start by using the np.diff function to compute the differences between consecutive elements. Then you can find where the sign of these differences changes (i.e. from positive to negative for maxima, and from negative to positive for minima). Then you can store the corresponding radius values for those indices and return them as the output of the function.

Task 6: Animating a 3D Wave Surface

For a challenge task, try to create an animation of the wave surface over time. However there is a slight twist: the equation for the wave surface now accounts for time decay as well, so the wave surface will not only change with time, but it will also decay over time. The new equation is: \[f(x, y, t) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r} \cdot e^{-\beta t}\] Where \(\beta\) is the decay constant of the wave with respect to time. You can use Matplotlib’s animation capabilities to create the animation. Here is a sample usage of the function:

>>>
x_mesh, y_mesh = create_meshgrid(10, 100)
>>>
animate_wave_surface(x_mesh, y_mesh, A=1.0, k=2.0,
...
omega=1.0, alpha=0.1, beta=0.05,
...
duration=10, fps=30)
Figure 1: Animated 3D Wave Surface Plot
Tip

Use the following incomplete code as a starting point for your animation function. You will need to complete the update_image function and the animate_wave_surface function to create the animation.

from matplotlib import animation
def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
    """Updates the wave surface for each frame of the animation.

    Parameters
    ----------
    frame : int
        The current frame number.
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    ax : matplotlib.axes._subplots.Axes3DSubplot
        The 3D axis object to update.
    """
    wave_surface = # write your code here!

    ax.clear()
    ax.plot_surface(x_mesh, y_mesh, wave_surface)
    ax.set_xlabel('X-axis')
    ax.set_ylabel('Y-axis')
    ax.set_zlabel('Wave Surface')
    ax.set_title('Animated 3D Wave Surface Plot')

    return [ax]



def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
                         k: float, omega: float, alpha: float, beta: float,
                         delta_t: float, num_steps: int):
    """Creates an animation of the wave surface over time.

    Parameters
    ----------
    x_mesh : np.ndarray
        The meshgrid of x coordinates.
    y_mesh : np.ndarray
        The meshgrid of y coordinates.
    A : float
        The amplitude of the wave.
    k : float
        The wave number.
    omega : float
        The angular frequency.
    alpha : float
        The decay constant of the wave with respect to distance.
    beta : float
        The decay constant of the wave with respect to time.
    delta_t : float
        The time step for each frame of the animation.
    num_steps : int
        The number of steps in the animation.
    """
    
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

    #| decorate the plot with labels and title here!

    anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
                                   fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
    plt.show()